##########################################################################
#####################        ## DESCRIPTIVES ##       ####################
#####################           Max Joosten         ######################
#####################   Last updated: 30.10.2023    ######################
##########################################################################

############################################################################
############################################################################
# Load packages
library(tidyverse)
library(sandwich)
library(stargazer)
library(sjPlot)
library(ggeffects)
library(writexl)
library(plm)
library(lme4)
require(AER) # for coeftest
library(vtable) # for summary table

windowsFonts(Times=windowsFont("Times New Roman"))

############################################################################
############################################################################
# Set working directory
setwd("")

############################################################################
############################################################################ 
# Load data
load("ref_v1.RData")
# Load data (only companies that have complete data)
load("ref_complete_v2.RData")

############################################################################
############################################################################
# TABLE 2: Distribution stringency
df %>% group_by(stringency) %>%
  summarise(n = n()) %>% t()

############################################################################
############################################################################
# TABLE 3: Descriptives
labs <- data.frame(name_df = c("id", "company", "sector", "cntry", "year",
                               
                               "overall_mean_addup_v2",
                               
                               "bmodel_mean_v2",
                               "i_policy_supplyselect", "i_policy_supplyend", "i_csr_exec_comp", "i_csr_team",
                               
                               "policy_mean_v2",
                               "i_policy_hr", "i_policy_hr_care", "i_policy_cl", "i_policy_fl", "i_policy_sc_hs_train", 
                               "i_policy_sc_esg_train", "i_policy_sc_hs", "i_policy_sc_hs_imp",
                               
                               "voluntary_mean",
                               "i_voluntary_ft", "i_voluntary_eti", "i_voluntary_eiti", "i_voluntary_gc",
                               "i_voluntary_asi", "i_voluntary_bci", "i_voluntary_fla", 
                               "i_voluntary_gcp", "i_voluntary_rspo",
                               
                               "sales_log", "intsales_log", "employees_log", "roa", "debttoassets"),
                   
                   label_df = c("Company ID", "Company", "Sector", "Country", "Year",
                                
                                "Overall (Mean)",
                                
                                "Business model (Mean)",
                                "Selection of suppliers based on HR/labour considerations", "Ending contract based on HR/labour considerations", "Executive compensation linked to CSR", "CSR Team",
                                
                                "Policy (Mean)",
                                "Human rights policy", "Human rights guaranteeing care policy", "Child labour policy", "Forced labour policy", "Supply chain health and safety training", 
                                "Supply chain ESG training", "Supply chain health and safety", "Policy Health and Safety Supply Chain Improvement",
                                
                                "Voluntary (Mean)",
                                "Member Fair Trade policy", "Member ETI", "Member EITI", "Member Global Compact", 
                                "Member ASI","Member BCI", "Member FLA", 
                                "Member GCP", "Member RSPO",
                                
                                "Sales (logged)", "International sales (logged)", "Employees (logged)", "Return on assets (t-1)", "Debt to assets (t-1)"))

st(df, vars = c("overall_mean_addup_v2",
                
                "policy_mean_v2", 
                "i_policy_hr", "i_policy_hr_care", "i_policy_cl", "i_policy_fl", "i_policy_sc_hs_train", "i_policy_sc_esg_train", "i_policy_sc_hs", "i_policy_sc_hs_imp",
                
                "bmodel_mean_v2",
                "i_policy_supplyselect", "i_policy_supplyend", "i_csr_exec_comp", "i_csr_team",
                
                "voluntary_mean",
                "i_voluntary_ft", "i_voluntary_eti", "i_voluntary_eiti", "i_voluntary_gc",
                "i_voluntary_asi", "i_voluntary_bci", "i_voluntary_fla", 
                "i_voluntary_gcp", "i_voluntary_rspo",
                
                "sales_log", "intsales_log", "employees_log", "roa", "debttoassets"),
   summ = list(c('mean(x)','sd(x)','min(x)','max(x)', 'notNA(x)')),
   summ.names = list(c("Mean", "SD", "Min.", "Max.", "N")),
   labels = labs,
   file = "../05. figures/sumstats")

############################################################################
############################################################################
# FIGURE 1: 
df_ay %>% ggplot(aes(x = year)) +
  geom_smooth(aes(y = overall_mean_addup_v2), linetype = 2, colour = "black", method = "loess", linewidth = 1.25) +
  theme_bw() +
  scale_x_continuous(breaks = c(2010:2020)) +
  labs(x = "", y = "Average RBC, by year") +
  theme(axis.text = element_text(size=20),
        axis.title = element_text(size=20),
        plot.title = element_text(size = 20, face = "bold"),
        text = element_text(family="Times", size=12),
        axis.text.x = element_text(angle = 90)) 

# Save FIGURE 1
ggsave("../05. figures/JoostenFig1.tiff", width = 12, height = 8, dpi = 300)

############################################################################
############################################################################
# TABLE 4
# Pooled 3FE OLS models (exact same as 'reg overall_mean_addup_v2 stringency employees_log roa_lag dta_lag voluntary_mean v2csreprss i.sector i.year i.cntry, vce(rob)' in STATA)
m1 <- lm(overall_mean_addup_v2   ~ stringency + employees_log + roa_lag + dta_lag + voluntary_mean + v2csreprss + as.factor(cntry) + as.factor(year) + as.factor(sector), data = df)
m2 <- lm(policy_mean_v2   ~ stringency + employees_log + roa_lag + dta_lag + voluntary_mean + v2csreprss + as.factor(cntry) + as.factor(year) + as.factor(sector), data = df)
m3 <- lm(bmodel_mean_v2   ~ stringency + employees_log + roa_lag + dta_lag + voluntary_mean + v2csreprss + as.factor(cntry) + as.factor(year) + as.factor(sector), data = df)

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1,
          m2,
          m3,
          se = list(r1,
                    r2,
                    r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05. figures/table_4.html",
          title = "TABLE 4. OLS 3FE results of Responsible Business Conduct activities at the firm level",
          omit = c("cntry", "year", "sector"),
          covariate.labels = c("Stringency", "Employees (log)", "Return on assets (lagged 1 year)", "Debt to assets (lagged 1 year)", "Voluntary score (mean)", "CS Freedom (Home state)"),
          dep.var.labels = c("Overall", "Policy", "Business Model"),
          add.lines=list(c("Year FE", "Yes", "Yes", "Yes"),
                         c("Sector FE", "Yes", "Yes", "Yes"),
                         c("Country FE", "Yes", "Yes", "Yes"),
                         c("Country-years", "339", "339", "339")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), c("df", "df_ay")))

############################################################################
############################################################################
# TABLE 5. Bottom / Middle / Top 
m1 <- lm(overall_mean_addup_v2   ~ stringency + employees_log + roa_lag + dta_lag + voluntary_mean + v2csreprss + as.factor(cntry) + as.factor(year) + as.factor(sector), data = df_ay)
m2 <- lm(overall_mean_addup_v2   ~ stringency + employees_log + roa_lag + dta_lag + voluntary_mean + v2csreprss + as.factor(cntry) + as.factor(year) + as.factor(sector), data = subset(df_ay, bot2010 == 1))
m3 <- lm(overall_mean_addup_v2   ~ stringency + employees_log + roa_lag + dta_lag + voluntary_mean + v2csreprss + as.factor(cntry) + as.factor(year) + as.factor(sector), data = subset(df_ay, mid2010 == 1))
m4 <- lm(overall_mean_addup_v2   ~ stringency + employees_log + roa_lag + dta_lag + voluntary_mean + v2csreprss + as.factor(cntry) + as.factor(year) + as.factor(sector), data = subset(df_ay, top2010 == 1))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))

# Regression
stargazer(m1,
          m2,
          m3,
          m4,
          se = list(r1,
                    r2,
                    r3,
                    r4),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05. figures/table_5.html",
          title = "TABLE 5. OLS 3FE results of Responsible Business Conduct activities at the firm level, bottom 20% / mid 60% / top 20% (in 2010)",
          omit = c("cntry", "year", "sector"),
          covariate.labels = c("Stringency", "Employees (log)", "Return on assets (lagged 1 year)", "Debt to assets (lagged 1 year)", "Voluntary score (mean)", "CS Freedom (Home state)"),
          dep.var.labels = c("All","Bottom 20%", "Middle 60%", "Top 20%"),
          add.lines=list(c("Sample", "Full", "Bottom 20%", "Middle 60%", "Top 20%"),
                         c("Year FE", "Yes", "Yes", "Yes", "Yes"),
                         c("Sector FE", "Yes", "Yes", "Yes", "Yes"),
                         c("Country FE", "Yes", "Yes", "Yes", "Yes"),
                         c("Country-years", "294", "294", "294", "294")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), c("df", "df_ay")))